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ABSTRACT 


Periodic solutions of the equation 

z + z + ( T - R + Rz 2 ) z + Tz = 0 

are studied numerically and, for large R and T , analyti- 
cally. The periodic solutions are unstable in a small strip 
of the (R,T) plane whose boundaries are determined. 


§ 1. Introduction 


In a recent paper (MOORE and SPIEGEL, 1966, henceforth 
referred to as MS) a thermo-mechanical oscillator was con- 
structed with the idea of studying, in a simplified system, 
the properties of non-linear overstability. A small buoyant 
element was allowed to move in a temperature-stratified fluid 
under the action of a linear restoring force. The buoyancy 
of the element was made to depend in a linear way on its 
temperature which in turn depended on the temperature of the 
surrounding fluid through Newton's law of cooling. The fact 
that a finite time is needed for temperature adjustment is 
essential to the functioning of the model. 

This cooling time provides a suitable time scale and 
if time is measured in these units and if a temperature 
stratification is chosen so that the surrounding fluid is 
unstable for |z| <1 and stable for jz| >1 (a cubic 

law is the simplest analytic function satisfying these re- 
quirements and this was used) the equation satisfied by the 
position z(t) of the small element is 

+ (T - R + Rz 2 ) — + Tz = 0 . (1.1) 

dt 3 dt 2 dt 
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R and T are large when the thermal dissipation involved 
in the cooling law is small. They can be thought of as the 
square of the ratio of the cooling time to the characteristic 
convective time in the absence of the spring and the square 
of the ratio of the cooling time to the period of free oscil- 
lations under the action of the spring alone. 

Since energy is available to the small element only for 
|z | <1 we would expect the motions to be bounded, whatever 

the initial conditions. Moreover, since the system is dis- 
sipative we would expect the system to enter a limit cycle 
for large values of t . In MS equation (1.1) was studied 
by time integration and it was found that, while a limit 
cycle was usually entered, there was a narrow strip of the 
( R, T) plane in which solutions were not asymptotically 
periodic. This strip is shown in Figure 1. The purpose 
of the present paper is to explain the existence of this 
strip. 

In MS it was suggested that the aperiodicity was due 
to the fact that when R and T lay in the aperiodicity 
strip the corresponding periodic solutions of (1.1) passed 
close to the singular point of (1.1) in its phase 
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(z, — , -- ^ ) space, enabling small numerical errors to 
dt dt 2 

push the phase point from one periodic solution to another. 
(Though the vacillation is induced by numerical errors, the 
same effect is produced on the actual physical system by 
small random disturbances. Thus the vacillation is a real 
physical effect. A rigid simple pendulum with just enough 
energy to make complete revolutions and subject to small 
random forces would display a similar vacillation.) Evidence 
was presented to support this theory of the origin of the 
vacillation and in particular a few periodic solutions of 
the required type were discovered using a relaxation proce- 
dure. 

In this paper a systematic search for the periodic so- 
lutions is made. In § 2. an analytic approximation to the 
periodic solutions of (1.1) is developed which is valid when 
R and T are very large. In this case the motion is ap- 
proximately adiabatic and satisfies the energy equation 

where 

6=1- T/R (1.3) 

and where E is the energy and B the buoyancy of the 
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element relative to the fluid at z = 0 — both are con- 

served in the absence of dissipation. To each E and B 
there corresponds a unique energy curve and to each such 
energy curve (apart from a single critical curve) corre- 
sponds a periodic solution. This degeneracy in the adia- 
batic theory can be removed in a familiar manner by con- 
sidering the average effect of the dissipation over the 
periodic solutions supplied by (1.2). This is valid since 
there are many oscillations per cooling time if R » 1. 

If we do this we find (§ 2a.) that 


= F(E,B,5) 

dt 

^ = G (E, B, 6) 

dt 


(1.4) 


where F and G are averages of the dissipative forces 
over the energy curves defined by (1.2). If F and G 
both vanish for particular values E 0 and B 0 then integra- 
tion of (1.2) with E 0 and B 0 inserted will yield an ap- 
proximation to the periodic solution of (1.1). This approxi- 
mate theory shows that a pair of distinct periodic solutions 
exists which pass vanishingly close to the origin in the 
phase space of (1.1) as 6 -» § from below. These solutions 
have a characteristic 'spiky' appearance with long flat 
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~A 


portions of small slope between the spikes. For 6 ^ 
they disappear. This is in good agreement with the time 
integrations which showed that, to very good approximation, 
the solutions were asymptotically periodic to the left of the 
line 6 ~ f in the (R,T) plane and had aperiodicity of 
the "jumping" type just to the right. 

Numerical solutions of (1.1) at very large values of 
R and T , such that (R,T) is in the aperiodicity strip 
have revealed that near the right hand boundary the aperi- 
odicity takes on a rather different form. Rather than a 
jumping between periodic solutions with different shapes 
and periods one has a continuous but aperiodic modulation of 
a single oscillation (see Figures 2 and 3) . Characteristic 
of the modulation is a gradual decrease in the amplitude of 
excursions on one side of z = 0 followed by a sudden jump 
to a nearly symmetric motion. The time between the sudden 
jumps is apparently random (at least we can determine no 
rules from inspection of about twenty time integrations in 
this region of the (R,T) plane) and it had no systematic 
dependence on R . 

An explanation of this effect in terms of the (E,B) 
plane is proposed in § 2d. The singular points of (1.4) 
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turn out to be all unstable if .75 > 6 > .62 and the 
point (E,B) traverses a limit cycle in the (E,B) plane. 

Thus instead of an oscillation of constant amplitude one 
gets a modulation of the amplitude (and shape) of the oscil- 
lation as (E,B) goes round its limit cycle, which is tra- 
versed in 0(1) cooling times. The sudden jump is associ- 
ated with the singular energy curve of the family (1.2) 
which separates energy curves of two types. It is conjec- 
tured that the aperiodicity of the modulation is also as- 
sociated with the singular curve. 

At finite R and T we can only hope to find the 
periodic solutions of (1.1) numerically. However, if R 
and T lie in the aperiodicity strip these periodic solu- 
tions are unstable and will not emerge from a time integration. 

In § 3a. we describe a relaxation method of finding the 
periodic solutions which works even when the solution being 
sought is unstable. We guess the shape of the solution 
and its period and improve the guesses by substituting them 
in the equation and calculating corrections. If we were pre- 
pared to solve non-linear algebraic equations for the values 
of the corrections at the grid points into which we divide 
the period (normalized to unity, so that the physical period 
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enters the equation as an eigenvalue) we could find the 
corrections in one step. Instead we assume formally that 
the corrections are small and solve the resulting linear 
equations. The process is repeated until the corrections 
fall below some pre-assigned value. 

The solutions found in this way agree well with the 
results of time integrations in the stable region of the 
(R, T) plane and with the asymptotic theory for R -* CO. 

In particular periodic solutions passing close to the origin 
of the phase space of (1.1) are found near the line 6 = 1 . 

In Floquet theory one studies the stability of peri- 
odic solutions of a non-linear ordinary differential equa- 
tion by finding out if small perturbations to the solution 
grow or decay. The similarity to our numerical procedure 
is obvious and our procedure yields the Floquet multipliers 
associated with the periodic solution as well as the peri- 
odic solution itself. It is found that the periodic solu- 
tions lose stability as a critical curve in the (R,T) 
plane is crossed from right to left and this curve is in 
good agreement with the right-hand boundary of the aperio- 
dicity strip. 
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§ 2 . The Method of Averaging 

2a. The averaged equations. 

In this section we take up the problem of finding ap- 
proximations to the periodic solutions of (1.1) in the case 
where R and T are large and R = 0(T) . 

We recall from the physical nature of the system de- 
scribed by (1.1) that infinite R and T correspond to 
no thermal dissipation, so that in this case the system 
will execute oscillations of constant amplitude. We shall 
call these adiabatic oscillations. If R and T are 
large but finite the dynamical time scale 0(R ) which 
characterizes the oscillations is very much shorter than the 
cooling time 0(1) which characterizes the dissipative 
process. Thus during one oscillation the dissipation can 
produce only a small effect and we can anticipate that the 
amplitude, degree of asymmetry and phase of the oscillations 
will change only very gradually. Consequently one can al- 
low for the effect of the dissipative terms by calculating 
their average over one oscillation (BOGOLIUBOV and MITRO- 
POLSKY 1961, p. 387 ff.) 

Let us introduce a new dimensionless time t = R 2 t , 
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so that T measures time in dynamical units. Then (1.1) 
can be written 


d 3 * 

cl T * 


(S-Z z ) 


dz 
d T 




( 2 . 1 ) 


The small dissipative terms are on the right and the adia- 


batic motions thus satisfy 


d 3 1 

d T 3 


(s-z 1 ) 



( 2 . 2 ) 


This equation was discussed in MS. If we integrate it twice 
we find that 


J. 

2 




+ ~ z 4 + = E 


(2.3) 


where E and B are constants of integration. Physically, 

E is the energy and B the buoyancy defect relative to 
the fluid at z = 0 of the oscillating element. The curves 
of constant E in the (z t ) phase plane are closed 
for all B , so that except for a critical curve which enters 
a saddle point, every curve, and hence every pair of values 
(E,B) , corresponds to a periodic solution of (2.2). 

Thus for any value of $ there is a doubly infinite 
family of adiabatic oscillations and it is up to the dis- 
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sipative process to decide which members of this family 
neither gain or lose energy and buoyancy over one cycle. 

These particular adiabatic oscillations will be our approxi- 
mation to the periodic solutions of (1.1) for the given 
value of 5 . 

We can integrate (2.3) again in terms of Jacobi ellip- 
tic functions but let us write the solution of (2.3) in the 
form 

z = z A ( t+0,E,B, 5) (2.4) 

where the new constant of integration 0 is a phase, and 
where the subscript is to emphasize that z a (t) is a 

solution of the adiabatic approximation (2.3) . 

We have suggested that for large R the effect of the 
dissipation will be to slowly change the shape of the adi- 
abatic oscillations. To allow for this we regard 0 , E 
and B not as constants but as varying only over many 
(actually 0(R ) ) cycles and we try to choose them so 
that (2.1) is satisfied in an average sense. As a conse- 
quence of the time dependence of the constants, 

d 2 - 3 3 « 324 c . 

dT “ 3T 3© u 3E C 35 
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where a dot denotes -i±- . Now if we insist that the last 

dT 

three terms vanish for all t we can ensure that 

21 - 3 Z a 
3T “ dT 


and by a similar restriction on the second derivatives we 
can ensure that 


JT ! 3 T l 

and so 

= i!ii „ ii ** e + ^ r + ^ p 

st j a T’ ar*30° st^c 11 8T*aB D 

Thus we have 


-* A 0 + =0 


de 


ae ae 


1^ + = 0 


> ( 2 . 5 ) 


h r t R = 

jT l ae° >t 1 5e 1 aT l ^6 u 
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where the final equation of the set is derived by insist- 
ing that z(t) satisfies the full equation (2.1). 

The second and third derivatives of z A in (2.5) can 
easily be expressed in terms of the first derivatives by 
using (2.3). For example, if we differentiate (2.3) with 
respect to E we find 


3 % 

3T3E ar 


I - 


111 

d E 






and using this and similar easily derived relations we can 
show that the solution of (2.5) is 




d*iiA 


(r + aT* 


)(** 


dl 


a + 32 a 


a E a 6 



d ?A 
dT 


*\ 


R y 'E = 


_ _ 2 I M 

? A B + Z A ~ 3 Z A 


> (2.6) 


w A -^z A 3 


J 


So far the analysis is exact and we could in principle 
solve these equations for 8 , E , B using the explicit 
solution (2.4). Needless to say, this is impossible and 
we replace the right-hand sides in (2.6) by their average 
values over one oscillation. Let us denote this average 
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by an overbar. Now 


d f l 


6 * A 



+ B =0 


and for a periodic solution 




so that, on taking the averages of the right-hand sides of 
(2.6), we can simplify the equations to 



These equations are considerably simpler in that there is 
no explicit time dependence and a solution is feasible. 

In particular we observe that since the average value of 
a power of z K is independent of 9 , the last two equa- 
tions suffice to determine E ( t) and B(t) . Thus the 
phase drops out of the problem (a reflection of the fact 
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that our system is unforced) and we are left with the pair of 
equations 


R l/l E = F(E,B,«) 
B = G(E, B, B) 


> ( 2 . 8 ) 

J 


where 

F(E,B,S) = + 


and 

G(e,b,*; - z a 


(2.9) 


( 2 . 10 ) 


We are interested in values E 0 and B Q which make 


F(E OJ M)=0 

and 

G(E o ,B o ,S) = 0 • 




r (2.ii) 


j 


We can solve (2.11) directly, but it is instructive to note 
that solutions of (2.11) are singular points of (2.8) in 
the (E,B) plane. This suggests that the structure of the 
solutions of (2.8) in the (E,B) plane will be of interest. 
For example, if the singular point (E 0 ,B 0 ) is unstable we 
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should not expect to find the corresponding adiabatic oscil- 
lation in a time integration of (1.1) . The idea of discuss- 
ing the nature of periodic solutions of a differential equa- 
tion in terms of the phase plane of the adiabatic constants 
is due to Andronow and Witt (STOKER 1950, p. 153) , who ex- 
amined second-order forced systems. 

Suppose that we are close to a singular point so that 


E- E/E 

B = 8,, » B 

then approximately 

~ _ (o _ foicr 

e * r E e + f b b 
B - G e M E ♦ G* B 




J ( 2 . 12 ) 


and the nature of the singular point is decided by the 
values of the four partial derivatives, unless the deter- 
minant 




(O) 

B 


-- (0) <o) 

' B 



(2.13) 
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(STOKER 1950, p. 44). When (2.13) is satisfied there is 
another solution of (2.8) near (E 0 ,B 0 ) so that we have a 
bifurcation. (The superscript means that the partial deriva- 
tives are evaluated at (E 0 ,B 0 ) ). 

The discussion of the structure of the (E,B) plane 
is greatly facilitated by the introduction of a generating 
function whose existence was pointed out to us by Whitham 
(1966 Private communication) . Let 


£) = § dr = f i„ii 


(2.14) 


where the integral is taken around the energy curve defined 
by E and B . Clearly H(E,B,5) is just the area under 
the energy curve. Whitham shows that F and G can be 
expressed in terms of H and its derivatives with respect 
to E and B and a short calculation gives 


F(E,8,6) 


P i 


1 ft dH / 3 -HB 
3S ~ V I -8 





and 

G(e,&,$) * " p ag f 


> (2.15) 

J 
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* . "l-> v» 

w nt; -L c 


(2. 16) 


p- M 
r d E 

is the period of the adiabatic oscillation. From (2.15) and 

( 2 . 11 ) 


I 2 H _ T- 

P 36 “ ft * . 


(2.17) 


The system of equations (2.11) has one trivial solution 
which we can dispose of before proceeding. We can verify 
analytically that 


I i fn litn -3W - o 

E-*0+ B-+0 d& 


I im 

E->0 


|iw AH. - no 

B"»0 3E 


while 


H(0,0) is finite. 


Thus E = 0 and B = 0 is a solution of (2.11) for all * . 

The corresponding energy curve is the critical figure of 

eight curve entering the saddle point (0,0) in the 

Cz. da ") plane, so that the only periodic solution with 
d t 
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E = 0 and B = 0 is z A = 0 , that is, no motion at alii 
This corresponds to the position of unstable equilibrium 
at z = 0 for the original equation (1.1). 

We next prove that (2.11) has non-trivial solutions 
only when 6 < # . Clearly H > 0 and, because the energy 

curves increase in size as E is increased, > o . Now 

9E 

an energy curve with E = 0 passes through the origin in 
the Cz, ^5- j plane. Thus an energy curve with E < 0 
is either entirely in the region z > 0 or entirely in the 
region z < 0 . In either case z A d 0 so that in view 
of (2.17) 

# 0 if E < 0 . (2.18) 

Thus solutions are possible only if E i. 0 and the condi- 
tion for a non-trivial solution to be possible follows at 
once from the form of F given by (2.15). 

We stress that we have not proved that periodic solu- 
tions of (1.1) are impossible at large R and T if 6 '> f 
only that, in this range, the periodic solutions are not 
approximated by periodic solutions of (2.3). This is con- 
firmed by the time integrations, since for 6 > ^ we found 
that the equation had periodic solutions, but that these 
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solutions had several maxima and minims per period, whereas 
periodic solutions of (2.3) have just one maximum and one 
minimum per period. 

h further transformation is useful. Let 


z = 5 s , 


E* = E/6 2 , (2.19) 

and 

B* = B/6 ^ . 

Then 

H = 2 VT J (E*-d"i S A ) V *is < 2 - 2 °) 

where s 0 and s, correspond to the points in which the 
energy curve cuts the z axis. We note that H depends on 
6 only through the factor 6 2 . Thus the curve T in the 
(E* , B*) plane on which =0 is independent of 5 , 

so that in view of (2.17) G and z" A vanish at every point 
of this curve. In general F will be non zero except at 
a finite number of points of T and as 6 changes these 
points, which are the singular points of (2.8) in the (E*,B*) 
plane will traverse F (Figure 7) . 

The symmetry of H with respect to B* shows that 
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-Jr. = 0 on B* = 0 (E* > 0) + (2.21) 

O-D 

so that one branch of T is the positive E* axis. Indeed 
it is obvious that z A vanishes on any symmetric energy 
curve — physically, the dissipation does not destroy the 
symmetry of an initially symmetric motion. These symmetric 
oscillations are simplest to discuss and we consider them 
in detail in § 2b. and § 2c. We return to the general case 
in § 2d., where we find that V has a branch on which 
B* -l 0 , giving rise to asymmetric oscillations. 

2b. The symmetric oscillations 

If B* = 0 the energy curves are symmetric about the 

axis. The critical curve now passes through the origin 
dT 

and corresponds to E* = 0 ; for E* > 0 the energy curves 
have a single symmetric branch (Figure 4a) . 

We now easily see that 

-So 

and 

If E* < 0 the energy curve is entirely above or below the 

r)H 

z axis and / 0. H(E*,B*) is not differentiable at (0,0) 
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f a 

I dH _ rz \ ds 

TWW* dE* v J^vW-shu *♦«,*) 

where 

C = 3 [I *(I*|E‘) V *] 

and 

S.* = 3 [(l*| E') l/2 - I ] . 

The integrals are easily evaluated (Byrd and Friedman 1954, 
p» 45) and we find that the equation 

/ir* £jj — 3 " jt !L u 
aE* ~ i - s H 

which will determine the value of E* corresponding to a 
steady oscillation is equivalent to 


"\ 

> ( 2 . 22 ) 
J 


HE* 

C i * % t') H 


l-HS 
I - S 



(2.23) 


where 


4 


i 




(2.24) 


and where, in the standard notation, E (k) and K(k) are 
complete elliptic integrals. A rather more convenient form 
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may be obtained by applying the imaginary modulus transfor- 
mation to E and K and after some manipulation we find 


/ A XO + 0 + 


(2.25) 


where 



Vi 1 / 

f ( | + sir»*ey*<ie 

O 

J ( i 


Clearly as uu increases from 0 to OO , $ increases mono- 

tonically from 1 to OO , so that as E* increases from 0 to 
OO , the right-hand side of (2.25) decreases monotonically 
from OO to 1 - 4/3 6 (we assume 0 < 5 < 3/4) whereas the 
left-hand side increases monotonically from 1 - 2/36 to 
OO . Thus there is a unique E* = E*(6) for which (2.25) 
is satisfied and the function E*(6) is easily found nu- 
merically. Once E* ( 6 ) is determined, s 0 and s 2 can 
be found in terms of 6 from their definitions (2.22) and 


the oscillation is given by 



and its period P bv 

— O v • 


1 iriA 





» 


P SYM is shown as a function of 6 in Figure 5 and we 
note that it increases rapidly as 5 approaches the criti- 
cal value 3/4. This is because as 5 -» 3/4 - , E*(5) -» 0 + , 

so that the energy curve is tending to merge with the figure- 
of-eight curve which has infinite period. In fact it can be 
shown by expanding for small E* that 



♦ 0(0 


as 6 -* 3/4 - . 


The form of the oscillations for 5 = 0.748 is shown 
in Figure 6. The flat portion of small slope which is de- 
veloped when 6 is near the critical value is due to the 
slow motion of the particle near the saddle point (0,0) 

Cz K , — J - ^ plane. This corresponds to a close approach 

' v dT / 

to the origin the phase space (z, ^ - 2 of the origi- 

v dt dt 2y 


nal equation (1.1) . 


The above analysis is restricted to the case R » T 
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for which 6 > 0 . If 6 < 0 very similar analysis can 
be given and it appears that there is again a unique sym- 
metric solution. 

The time integrations showed that there was a stable 
symmetric solution only when 6 < . 54 , approximately. This 
suggests that the symmetric solution is unstable when 

.54 < 6 < .75 

This is confirmed in § 2c. The time integrations showed a 
stable asymmetric oscillation when 5 was slightly greater 
than .54. We consider this in § 2d. 

2c. The stability of the symmetric solution. 

We have seen that a symmetric solution exists in the 
range - oo < § < 3/4. In this section we discuss its sta- 
bility by examining the nature of the singular point (E* f O) 
where E*(5) is the solution of the transcendental equa- 
tion (2.25) . 

It follows from the symmetry of F with respect to B 
th3t F e (E , 0/ 6 ) = 0 , (2.26) 

so that (2.12) assume the form 
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R = 


(®)*v 
w ' * p 

*" E “ 


B* = G e <0, E + G B (rt B 


and it follows from STOKER 1950, p. 44, that the point 
(E o ,0) is 


I) a stable node if 


f/ 0) g/ 0> > 0 ; F p - (<,) + G b (<>) < 0 


II) an unstable node if 


f/°'g b (<>> > 0 ; F E (6) + G 9 (o) > 0 


III) a saddle point if 


f e (0) g b (<>) < 0 


Now it can be shown from the results of § 2b. that 


F F (E (6) ,0,6) < 0 


(2.27) 


Thus case II cannot occur and we have a stable node if 

< 0 and a saddle point of > 0 . If G B *°* =■ 0 

there will be a bifurcation. 

We cannot calculate G B without considering asymmetric 
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solutions. However, since the departure from symmetry is 
small, we can use a perturbation method. Thus to calculate 
the function G g (E o ,0,6) we write 

z = z a (t + e,E o (6),0,6) + y . (2.28) 

and substitute in (2.3) , retaining only terms linear in y . 
Thus we find 

y - y(z A /z A ) = - B(z a /z a ) 

whose solution is 

y = - Bz a 

where D is a constant of integration. However the value 
of D does not alter the average values involved in the 
calculation of g b and we set it equal to zero in subsequent 
work. 

Substituting the explicit form of z A we find, after 
some algebra, that 


(z A /z ) dt + Dz. 


y = 


* w*v>’ 


{ 




(2.29) 


It is now a straightforward matter to calculate y and since 
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uiuci >-'j. uui 3 uXiiita Lluii j 


BG, 


«>) 



we find that 


r (0> * ij s «* U usj 1 - 

RS (s. l *s t 1 ) 2 [s. 1 “ 9 -/** 


Clearly the sign of G,,* 0 ' is decided by the sign of the 
quantity in the square bracket in (2.30) and substituting 
for s 0 its definition (2.22) we have 



if 


E'us) 


> i a 


From the numerical solution of the transcendental equation 
(2.25) which defines E*(5) we find 


G t M (£.,o,s ) >< o if {>, ,su- , 

and we conclude that the symmetric solution is stable only 
if 

6 < . 561 • • • 


.30) 
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This is in good agreement with the time integrations of 
(1.1) for large R . 

The critical value of 6 corresponds to a bifurcation. 
In the next section we discuss the structure of the (E*,B*) 
plane as a function of 6 and we show that a stable asym- 
metric solution bifurcates from the symmetric solution at 
the point corresponding to the critical value. This is also 
in agreement with the results of the time integrations. 

2d. The structure of the (E*,B*) plane 

We have seen that (2.11) have a symmetric solution for 
6 < 3/4 . To complete our discussion we must examine the 
possibility of solutions with B* ^ 0 . In view of the 
symmetry of the (E*,B*) plane with respect to B* we 
need consider only positive B* . 

The energy curves (Figure 4) are 



and we must first discuss how these energy curves depend 
on E* and B* . The positions of equilibria satisfy the 


cubic 



3s + 3B* 


0 


(2.32) 


s 3 - 


so that there are three positions of equilibrium if 
O < B* < 2/3 and one if B* > 2/3 . Moreover in the first 
case the central position is unstable and the outer posi- 
tions stable whereas in the second case the sole equilib- 
rium position is stable (these statements are proved in MS) . 
Thus if B* > 2/3 the energy curves are closed ovals, which 
increase in size as E* increases from the value E*(s e ) 
appropriate to the position of equilibrium s e . If B* < 2/3 
the situation is more complicated. Let s e ,0 ( s^ 2 ', s e (3) be the 
positions of equilibrium where 


10 


f*> 


( 3 ) 




s e < s e < S e 


so that s 

<0 

e 

and s e ' 3) are the stable positions. 

When 

B* = 0 s e 

(i) _ 

- >/3 , s e <2) = 0 and s e f3) = + V3 . 

As 

increases 

s 0) 

*e 

and Sg^ 1 decrease and s e * 2 * increases 

until when 

B* 

= 2/3 s e <*> and s e * 3 ^ merge at the 

value 

1 while 

s (,) 

b e 

- -2 . In Figure 7 we sketch the curves 

E* (s e <0 ) 

(and 

its continuation E*(s e ) ), marked 

CD , 


E*(s e (i ^) , marked CA , and E*(s e ***) marked OA . Suppose 

we fix B* < 2/3 and increase E* from E*(s e <l) ) so that 


we move along a line parallel to the E* axis. The energy curve 
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consists first of a single oval surrounding which 

increases in size as E* increases. After we cross CA 
the energy curve has two branches, one surrounding s e (l * 
and one surrounding s e (3 ^ . As E* increases both ovals 
grow until on the curve OA they merge and enter the un- 
stable point s e . The energy curve is now a single fig- 
ure-of-eight curve and the motions are non-periodic. This 
curve is called the separatrix. When E* increases still 
further, the energy curve is a single closed curve surround- 
ing the separatrix. 

In short, we have a double branched energy curve in 
the crescent shaped region C A O and a single branched 
energy curve everywhere else to the right of CD . The arc 
OA corresponds to the separatrix. 

When the energy curve is double branched, the branch 
surrounding s e is entirely in the region s > 0 so 
that clearly z" A ^ 0 on this branch. Thus when the curve 
is double branched the branch surrounding s e **' is the 
branch on which H is calculated. 

The next step in the discussion is to determine the 
curve r on which BH/&B* = 0 . We have already seen that 
F has a branch B* = 0 corresponding to symmetric oscil- 
lations. When B* > 0 solutions are sought numerically. 
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of r bifurcates from B* = 0 at the point Q (•which was 
already located in § 2c.) and returns to the symmetric branch 
at 0 , crossing CA at the point P . The dotted portion 
of F is only inferred since, owing to the proximity of T 
to the critical curve OA corresponding to the separatrix, 
the numerical work becomes very difficult. Fortunately, on 
the curve CA , BH/BB* can be expressed in terms of ele- 
mentary functions and the point P located exactly, provid- 
ing a check on the numerical calculations. The closeness of 
r to OA is shown by the fact that A has coordinates 
(1/4, 2/3) whereas P has coordinates (.24, .6572 ...) 

Near 0 an approximate analytic treatment becomes pos- 
sible and rather lengthy calculations (Appendix A) show that 
in the sector AOB* 


an _ 
a 6* 


Hsft8 Vl 



( B**- 2 £ 




(2.33) 


whereas OA itself is 

B * 3 = 2E* 


(2.34) 


so that BH/bB* vanishes on the curve 



The proximity of the curves F and 0 A is clear from this 
formula . 


Once r is determined the point on F at which 


4e 


* 9H 
BE* 


I - <? 


H 


can be found. The arrows on Figure 7 show the direction $ 
increasing. The bifurcation point Q is, as already seen 
in § 2d., 6 = .561 , the point P is § = 5/8 and at O , 
fi = 3/4 . The period P ASM of the asymmetric oscillations 
increases rapidly as 5 increases (Figure 5) ; the approxi- 
mate solution given in Appendix A shows that 


P ASM ~ irV?/[ 2(3-46)] 


(2.36) 


in contrast to the logarithmic behavior of P &y . M . 

For 6 just smaller than 3/4 the energy curve is just 
inside the lower loop of the separatrix. Consequently the 
phase point spends most of each period near the point s e ^ , 
traversing the rest of the energy curve in a time of 0(1) . 

As a result the oscillations have a spiked appearance, with 



long flat portions between the spikes. Thus the asymmetric 
solutions also correspond to a close approach to the point 
(0,0,0) in the (z,z,z) phase space of (1.1) when 5 -* 3/4 

The nature of the singular point (E*,B*) was also de- 
termined. On E*Q , the singular point is a stable node and 
on QO , it is a saddle point, so that stable symmetric oscil 
lations are possible only if 6 < .561 . On QPO , it is at 
first a stable node, then a stable spiral point and then at 
R , where 5 = .62 , it becomes an unstable spiral. Thus 
for .561 < 6 < .62 stable asymmetric solutions are pos- 
sible. In Figure 8a we show the result of an actual time 
integration when 5 =0.6 . The point (E*,B*) is spiral- 
ling into a stable spiral point as it should. 

When 5 > .62 all the singular points of the system 
(2.8) are unstable. This means that there are no possible 
stable oscillations in this case. If we integrate (2.8) 
in time E* and B* will not tend to a fixed point as 
t - 00 and it is plausible to suppose that there will be 
a limit cycle. This limit cycle cannot surround the saddle 
point on OQ , so it must surround either the unstable 
spiral on OR or the point 0 itself, which, as we saw 
in § 2 a. is always a singular point of ( 2 . 8 ). 
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If (E*,B*) is describing a limit cycle in the (E*,B*) 
plane the oscillation corresponding to E* and B* will 
undergo a periodic modulation of amplitude and shape. The 
period will be 0(1) cooling times and will depend only 
on 6 . Such periodic modulations of a fundamental peri- 
odic solution have been observed for forced motion of a 
system of two degrees of freedom (Hayashi 1964, p. 309). 

What is novel in our system is 

a) the sudden jump in amplitude 

b) the irregular nature of the modulation. 

The sudden jump is reasonably simple to explain. Sup- 
pose that when the limit cycle is surrounding the unstable 
spiral point on the asymmetric branch of T (that is, a point 
below R ) its lower portion crosses OA . If its lower 
portion crosses OA it is easy to prove that it does so 
from left to right. Now crossing OA in this direction means 
going from the crescent shaped region in which the energy 
curve has two distinct branches to the region where it sur- 
rounds the separatrix. Thus as (E*,B*) crosses OA on 
its limit cycle the corresponding energy curve jumps from 
an oval just inside the larger portion of the separatrix to 
a dimpled curve surrounding the whole separatrix (Figure 4b) . 
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Thus the amplitude and shape of the oscillation will change 
in just the manner observed in the time integrations. In 
particular the phenomenon will not occur for 6 < .62 , 
since the asymmetric solution is then stable, as we have 
already seen. This is in good agreement with the results 
of the time integrations. 

We have not been able to prove that the limit cycle 
crosses OA . Instead we have verified the conjecture by 
time integration of (2.8). The results for 5 = .625 are 
shown in Figure 8 b. The solution curve appears to be ap- 
proaching a limit cycle of the conjectured form. 

For 6 > .7 the limit cycle still crosses OA , but 
surrounds the point 0 , being symmetric about the E* axis 
This presumably reflects a change in character of the singu- 
lar point as it moves along F towards 0 with increasing 
6 . However we have not been able to compute r in this 
region, so this remains a conjecture. 

Finally, we consider the irregular nature of the oscil- 
lation. Suppose that a solution curve of (2.8) crosses OA 
On that part of the solution curve near the crossing point 
the energy curve corresponding to (E*,B*) is very close to 
the separatrix. Thus small changes in E* and B* can 

make a large change in H (recall that and are 

9E* SB* 
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infinite on OA ) . Thus (2.8) is very sensitive to errors 
in E* and B* and solution curves which start at neigh- 
boring points on the left of QA can widely diverge to the 
right of OA . Suppose we start numerical integration at a 
point on the limit cycle. After crossing OA the solution 
curve will be pushed off the limit cycle by the magnified 
errors. Subsequently it will spiral back toward the limit 
cycle. Thus it will again encounter OA and suffer a sec- 
ond random deflection. This process will repeat itself in- 
definitely so that the solution curve will never coincide 
with the limit cycle. 

This does not conflict with the capture property of 
stable limit cycles since this depends essentially on the 
fact that there is a unique solution curve through every 
point of phase space. The magnified errors can be thought 
of as random external forcing functions acting on the system 
(2.8) in a narrow region surrounding OA and in this region 
there can be many solution curves through any point. 

The irregular nature of the solution curve was verified 
by long time integrations of (2.8). 

This explanation of the aperiodicity of the modulation 
depends on the existence of an unstable spiral singular point 
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of (2.8). As we have seen the asymmetric singular point 
goes unstable at 6 = .62. This is in good agreement with 
the right-hand boundary of the aperiodicity strip found by 
time integration. 


§ 3. Periodic Solutions and the Aperiodic Range 

3a. The periodic solutions for finite R and T. 

In the preceding sections it was seen that for very 
large values of R and T there can exist symmetric periodic 
solutions of (1.1) so long as 5 < .75* These symmetric 
solutions are stable only in the range 6 < .561 and above 
this critical value an asymmetric solution exists, which 
itself is stable only in the range .561 < 6 < .62. This 
asymptotic theory thus predicts that in the range .62 < 6 < .75 
there will be both symmetric and asymmetric solutions, but 
that neither type will be stable. 

It is to be expected that the asymptotic theory will 
be approximately valid at large but finite R, and the above 
description accounts qualitatively for a number of the pro- 
perties of the solutions at R = 100 found by the numerical 
time integrations in MS. For example, it was found that 
the symmetric oscillation always occurred for small values 
of 6, but that the system settled into an asymmetric os- 
cillation when T was made less than about 44(6 > .56). 
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When T became less than 40(5 > .60) the oscillation became 
aperiodic. No periodic solutions were found in the range 
40 > T > 25 (.60 < 6 < . 75 ) and the periodic solutions 
which appeared for 6 .75 were of the multiple-peaked, 

or squegging, sort which cannot be predicted by the asymptotic 
theory since this type of solution is not exhibited by 
equation (2.3). Of course it is not possible to be sure, 
when a stable periodic solution is found, that it is the 
only one that exists; and unstable periodic solutions of 
the sort discussed in the preceding sections cannot be 
found by the straightforward time-integration technique. 

In this section we present examples of the symmetric 
and asymmetric solutions at finite R in their respective 
unstable ranges. By using special numerical techniques, 
it is possible not only to exhibit the unstable periodic 
solutions, but also to show by the application of Floquet 
theory that they are indeed unstable in an infinitesimal 
sense. This offers further confirmation that the "aperiodic" 
behavior discovered in MS is indeed vacillation among 
several truly unstable periodic solutions, and that the 
predictions of the asymptotic theory are qualitatively 
valid at large but finite R. 

Our numerical procedure depends upon a relaxation 
technique, in which it is demanded that the solution be a 
periodic one. The problem is changed from a one-point 
boundary (initial value) problem to a two-point boundary 
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problem with periodic boundary conditions. Thus only 
periodic solutions can be found (and of course there is 
no certainty that one has found all of these). The period 
appears in the role of an eigenvalue. 

In order to apply this procedure equation (1.1) was 
transformed into a system of three first-order equations, 
and the independent variable was taken to be x = tA, where 
A is the (unknown) period. The system thus becomes (with 
y x = z) 


dy 3 

dx 


y2 


dy2 

5T = y 3 


^3 

dx 


-[Xy 3 + (T-R+By 1 2 )X 2 y 2 + TxV, 


(3.1) 


with the boundary conditions 

y ± (o) = y.(i) i = 1 , 2 , 3 . (3.2) 

For definiteness, it is also necessary to fix the phase at 
some point; for example, we normally take 

Yl(0) = 0. (3-3) 

By a numerical procedure described in Appendix B we proceed 
to solve the system (3*1) subject to the conditions (3*2) 
and (3*3). One begins with a set of trial "solutions" 
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y.j_(x) and a trial value of 7 \, which may either "be guessed 
or taken from the solution of some previous case. In a 
series of iterations the functions and the eigenvalue are 
systematically corrected until the equations and conditions 
are fulfilled to high accuracy. The nature of the initial 
trial functions normally determines which solution, if any, 
is finally obtained. With experience it becomes relatively 
easy to obtain quite rapid convergence to the desired solution. 
The number of iterations required varies considerably, 
depending upon both the nature of the solution and the 
goodness of the trial solution. Regular, nearly sinusoidal 
solutions are very easily obtained, while strongly asymmetric 
solutions and those with sharp peaks require many iterations 
with gentle corrections. 

In order that the solutions of the difference equations 
may approximate as closely as possible those of the corre- 
sponding differential equations, the number of mesh points 
in the interval [0,1] of the independent variable was in- 
creased until it was determined that the results no longer 
changed when the number was further increased. It was found 
that 1000 points is usually a sufficient number. When this 
point is reached, the error in the solution presumably 
depends chiefly upon the accumulated round-off error in 
the arithmetical operations. In order to keep this as small 
as possible, double-precision floating-point arithmetic 
was used throughout (i.e., round-off of an individual number 
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occurs in the sixteenth significant decimal digit). It 
ues possible to use as many as 4000 equally spaced intervals, 
and in a few difficult cases it proved necessary to go to 
this number in order to achieve reliable results. 

Some of the results obtained with this procedure are 
reproduced in Figure 9* for R = 100 and various values of 
5 in the range -55 £ 6 £ .74* Time integrations for all 
cases of this figure are shown in MS. Figure 9a shows the 
single symmetric solution at 5 = .55* just before the bifur- 
cation point. At the slightly larger value 6 = .60 (Figure 
9b) the asymmetric solution has appeared and it is the stable 
one. Both this asymmetric solution and the symmetric one of 
Figure 9a agree with those found by time integration in MS. 

As 6 groxtfs the asymmetry becomes more pronounced as seen in 
Figure 9c, which is the solution at 6 = .61, Just into the 
range where the asymmetric solution has also become unstable 
and aperiodicity appears in the time integrations. Continuing 
through the unstable range (Figures 9d and 9e)* it may be 
noted that as in the asymptotic theory the flat "wings" 
in the symmetric solution become more pronounced, while the 
asymmetric solution gradually turns into a sharp peak fol- 
lowed by a long "stillstand". The period of the latter solu- 
tion increases rapidly, as predicted by the asymptotic 
theory. (All the asymmetric oscillations shown here are 
for B <0; of course the solutions for B >0 exist as well.) 
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At 6 = .74 (T = 26) we are nearing the end of the range 
where solutions exist according to the asymptotic theory, 
and the predictions of that theory are becoming less reliable. 
The symmetric solution (Figure 9f) is still very similar to 
that of the asymptotic theory (§ 2b). We were not able to 
obtain the asymmetric solution for 6 > .71; as the strongly- 
peaked character of this solution apparent in Figure 9© 
becomes still more pronounced, the first and second deri- 
vatives of the solution function become very large and the 
numerical computations come to be dominated by the round- 
off error. Thus we cannot be sure that the asymmetric solu- 
tion still exists at 6 = .74 for R = 100. There has also 
appeared a new type of periodic solution, the double-peaked 
oscillation shown at the bottom in Figure 9f • This seems 
to be an extreme case of the squegging behavior seen in 
MS at still lower values of T. We did not succeed in finding 
other periodic solutions at 5 = .74, such as a 3-peaked os- 
cillation or the asymmetric 2-peaked solution found at 
T = 25 (see the lowest curve in Figure 3 of MS), but it is 
quite possible that they exist. In any case it appears that 
there are no stable periodic solutions at 5 = .74, as shown 
by the extensive time integrations previously reported (see 
Figure 4 of MS). 

We succeeded in finding single-peaked symmetric os- 
cillations up to 6 = .748. For larger values of 6 no single- 
peaked periodic solutions were ever found, in good agreement 
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with our expectations on the basis of the asymptotic theory. 
The single-peaked symmetric solution at 6 = .748 and R = 100 
is shown in Figure 6, contrasted with that calculated from 
the asymptotic theory. The solution at R = 100 has a period 
about 5$ greater than that of the corresponding asymptotic 
solution and its shape is also somewhat different. The 
maximum occurs less than one-quarter period after the zero 
and there is a tendency toward a "stillstand" on the falling 
part of the curve. The nature of the solutions at 6 = .748 
is discussed in more detail below (§ 3o.) 

The periods of the symmetric and asymmetric solutions 
at R = 100 are compared in Figure 5 with those computed from 
the asymptotic theory. A study of this figure reveals the 
following points. The period of the symmetric solution agrees 
very well with that predicted throughout the entire range, 
differing slightly only when 6 becomes very close to . 75 , 
as noted in the preceding paragraph. The point of bifurca- 
tion occurs at a very slightly larger value of 6 (6 .566 

as compared with 6 . 561 in the asymptotic theory). Just 

above the bifurcation the period of the asymmetric solution 
at R = 100 is slightly lower than the asymptotic value, 
but it has become larger by the time 6 = .625 (the numerical 
asymptotic calculations were not carried beyond this point). 

In Figure 5 the period P ASM given by the approximate formula 
(2.36) is also shown. The periods at R = 100 are consistently 
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higher than P In the range .62 ^ 8 <: .70, the period 

of the asymmetric solution at R = 100 may be closely ap- 
proximated by the empirical formula 1/P = . 495(*738 - 6) 
as compared with (2.36) which gives 1/P = .468(.750 - 5). 

nujyi 

At 5 = .71, the highest value for which we were able to 
compute the asymmetric solution, the period appears to be 
increasing even more rapidly, possibly indicating that the 
asymmetric solution disappears at a value of 5 somewhat 
below .75. However the solution at 6 = .71 was obtained 
only with some difficulty and the computed period may not 
be reliable. 

The double-peaked solution on the right in Figure 9f 
was found in the range .73 < 5 < -75 • The periods of this 
solution are also indicated in Figure 5* 

In addition to the periods, the quantities E* and B* 
were evaluated for the periodic solutions at R = 100, allowing 
a comparison with the r curves of the asymptotic solution 
shown in Figure 7. For the symmetric solutions the dependence 
of E* on 6 given by Equation (2.25) was very closely repro- 
duced by the solutions at R = 100. Furthermore the segment 
QP of the F curve for the asymmetric branch agrees very well 
with that found at R = 100. (For reasons indicated in § 2d 
i\re do not have numerical computations of the asymptotic 
solution along P0; the portion of this segment near 0, where 
the approximate analytic treatment of Appendix A is applicable, 
is not accessible to our periodic numerical computations, as 
discussed earlier in this section. ) 


- 44 - 



3b. Stability of the periodic solutions. 

There is yet another point at which the numerical 


calculations can he compared with the theory; namely, we 
can investigate numerically the stability of the periodic 
solutions that have been found. This is accomplished by 
the application of Floquet theory (STOKER 1950, p. 193 ff . ) 

The problem investigated in this theory is whether, 
a periodic solution of a non-linear equation having been 
found, a small perturbation of this solution will remain 
bounded. Assume that we have found a periodic solution 

x ) 

xj of the system (3*1) with period "K. Consider 
x) 

now a vector differing only slightly from the above solution: 



y, + w “ 

y£ + wi where Iw^ (x) | « |y. (x) | . Then it is found 

_ y 3 + » 3 J 

2 

that upon neglecting terms of order w^ and higher, the 
v/. must satisfy a set of 3 linear equations 


dWj 

dx~ = f ±3 W J’ ^ i,J = 1,2,3 (3-*0 


with periodic coefficients: f. (x + X) = f. . (x). 

ij ij 

It is possible to find a noxmal fundamental set of 
solutions w 1 ? (x) of (3.4), and all possible solutions of this 
linear system can be expressed as linear combinations of the 

n ryr 

wj (x). The solution y 0 of (3-1) is said to be stable if 

LyiJ 

all solutions w j( x ) of (3*^) remain bounded as x ® which 
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will be true if and only if all the normal solutions are 
bounded. 

Now it is shorn in Floquet theory that all the normal 

solutions w* (x) can be expressed in the form 
1 ' 

n / \ a.? x n / \ 

w ± (x) = e i ©i (x) 


where the cp? (x) are periodic functions of x with period X 
and ou is given by 



The (in general complex) quantities are the Floquet 
multipliers and are the eigenvalues of the Floquet matrix 
A. This matrix is easily evaluated for any given periodic 
solution; indeed it is closely related to certain matrices 
which occur naturally in our iteration scheme to find the 
periodic solution and, once the scheme has converged, the 
matrix A is available with little further computation. The 
way in which A is obtained is discussed in detail in 
Appendix B. 

r y i ( x n 

It is now clear that the solution y„ (x) is stable 

Ly 3 (x)_ 

only if the real parts of all the are non-positive, for 
in that case all solutions of (3*4) will be bounded. Hence, 
we must evaluate the eigenvalues cx^ of A, and we shall expect 
the corresponding periodic solution to be stable only if 
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I I s 1 for all cr_^. Since it can be shown that one of the 
eigenvalues must always be equal to unity, it is the other 
two which must be determined. In practice, we solve numerically 
the cubic characteristic equation | ^ - alj | = 0, and verify 
that one of the roots is always equal to unity as a crude 
check on our procedure. 

In Figure 10 the quantity a = X” 1 in I a I, where 

1 max 1 

CT max is the Floc l uet multiplier having the largest absolute 
value, is plotted as a function of 6 for the symmetric and 
asymmetric solutions at R = 100. It is seen that the 
symmetric solution is stable (a < 0) for values of 6 below 
the point of bifurcation. For 6 £ .556, the symmetric 
solution is unstable and the asymmetric solution is now the 
stable one. This solution remains stable up to 6 « .606, 
after which it too becomes unstable. This is in excellent 
agreement with the results of the time integrations at R = 100 
and also agrees satisfactorily with the predictions of the 
asymptotic theory. 

In part of the stable range of the asymmetric solution, 

.589 ^ 6 £ .604, the Floquet multipliers are complex, and in 
this region a is always equal to -O.5. This region is indicated 
by the dashed part of the curve in Figure 10. The curve corre- 
sponding to the asymmetric solution in Figure 10 extends only 
to 6 css .64, although the solution was followed to 6 = . 71. 
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This is because the Floquet matrix becomes very ill-conditioned 
at the larger values of 6, and its eigenvalues can no longer 
be reliably determined numerically, even though the solution 
found fulfils the difference equations to a sufficient 
degree of precision. 

The quantity a for the double-peaked solution (Figure 9^) 
whose period is plotted in Figure 5 is also shown in Figure 10. 
This solution is never stable throughout the range of its 
existence, although it is similar, but not identical, to 
another double-peaked solution which is stable in a small 
range (see §3c). Finally, it will be remarked that the 
symmetric solution again appears to become stable for 6 > .745. 
This unexpected behavior is also discussed in the following 
section. 
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3c - Some properties of the solutions near 5 — 3/4. 

It was noted from Figure 10 that the single-peaked 
symmetric solution at H = 100 becomes stable when 6 s .7455* 
This was verified by a time integration at 6 = .748 > a portion 
of which is shown in Figure 11a. However this is not the 
only stable solution; with different initial conditions the 
apparently stable double-peaked solution of Figure lib was 
also found. This is the periodic solution found by time 
integration at 6 = .750 and reproduced in Figure 3 of MS. 

The solution curve is asymmetric, the peaks above the axis 
being higher than those below. There is also a tendency 
toward a "stillstand" before crossing the axis from above, 
but not from below. 

The symmetric double-peaked solution of Figure 9f also 
exists at 6 = .748 and is shown in Figure 11c. Its period 
is slightly less than that of the solution in Figure lib, 
and it is, according to the Floquet multipliers, unstable. 

It was verified that when a time integration is begun with 
initial conditions corresponding to the solution of Figure 11c, 
the oscillation starts off looking like this solution and 
then settles gradually into the stable solution of Figure lib. 
We did not succeed in finding the asymmetric double-peaked 
solution by the numerical method of § 9» bacause of the 
proximity of the more easily found symmetric double-peaked 
solution. 
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APPENDIX A i Approximate Solution near 0. 

Let s 0 , s 1# s s , s 3 be the roots in ascending order 
of the quartic 

s 4 - s 2 + B*s - E* = 0 (A . 1) 


These roots are the points in which the energy curve cuts 
the s axis and in the crescent shaped region QAC of 
Figure 7, they are all real. The branch which intersects 
the s axis at s x and s 2 has s > 0 at every point, so 
that if we are seeking solutions which make s = 0 , we must 
take our averages on the branch which intersects the s 
axis at s 0 and s T . Thus 



Standard transformations show that 


3 /* r * 

3H _ 4^/6 $ ; *S 6 f (i -«<,* sn*a) . 

>6 ' ycsj-s.xsri.) J 0 < ' -•<* ' 
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where 


c( 


z 





z 

I 


ii 

s„ 


t 


and 


IS} -s,)(s z -s 0 ) . 


(A. 4) 


The standard notation for elliptic functions and integrals 
is used throughout. 


When E* and B* are small we can solve (A.l) approxi- 
mately to find 


s 0 = - ^6* + B* 

s 1 = B* - Jb * 2 - 2E* 
s s = B* + >/B* 2 - 2E* 
S 3 >fs + B* 


(A. 5) 


The roots Sj and s 2 coincide when B* 2 = 2E* , so this 
is the approximate equation of OA when E* and B* are 
small. Our approximation is valid in the region AOB* so 
that B* 2 > 2E* and the roots are all real. 
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We can now easily find approximate expressions for «(* 
** and ft* , using (A. 5). In particular 

ft 1 - i-^ v /6 ‘ , - 2E ’ 

so that, since B* and E* are small the complementary 
modulus -ft , defined by 

■k * y/l - ft 1 

is small. In this circumstance the integral in (A. 3) can be 
evaluated approximately (Byrd and Friedman 1954, p. 301) and 
we find 


r K * a 

/ -«<, Sl> U 

I -•( i *n*u 
O 


du 




Substitution of the approximate expressions for •( * «( * and 
ft 1 * gives equation (2.33) of the text. The integrals for H 
and 3H/6E* are rather easier to estimate and we will not 
give the details of the algebra leading to (2.36). 
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Appendix B. A numerical method to obtain periodic solutions 
of nonlinear ordinary differential equations. 

In this Appendix we describe in detail the numerical 
procedure used to obtain the periodic solutions discussed 
in § 3* The method is based on one which is widely used in 
stellar structure calculations (L. G. Henyey, J. E. Forbes, 
and N. L. Gould, Astrophys. J. 139, 306 [1964]) and has been 
put into the form described here with the help of Dr. L. B. Lucy. 

Consider the Nth order system 



We are seeking periodic solutions of this system, namely a 
set of functions y (t) having the property 

y* (i ♦ A) = *ji (i) (b.2) 

for some real ~K and all t. As it is inconvenient in numerical 
work to have a variable interval for the independent variable, 
we change to a new independent variable x = t /~K and then 
confine our attention to a single period: 0 s x £ 1. Our 

system is now 

" h ( 1} ; M j * '< l ’ / A/ (b.3) 

and the required property (B.2) is expressed as a set of 
boundary conditions: 
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(B.4) 


— • , N. 

The problem now has the character of a two-point boundary 
problem, with the period \ as an eigenvalue. The condi- 
tions (B.4), however, do not suffice to determine a solution 
of (B.3), since the phase is not fixed; i.e., any solution 
can be arbitrarily translated in x and remains a solution. 
Thus we need some (quite arbitrary) condition to define the 
phase of one particular y.^. Normally we take 


U, (0) « C * corut. 


(B.5) 


It is important only that the value of the constant lie 

within the (normally bounded) range assumed by y-j_. 

In order to construct difference equations, we divide 

the range [0,1] into M intervals, not necessarily of equal 

length. Thus there will be M + 1 mesh points, labeled x , 

x 2 , x^, x^ +1 , with x = 0 and x^ + ^ = 1. The interval 

between two mesh points is then Ax 111 = x^-^ - Values 

of the y^, etc. at x = will be denoted by a superscript, 

thus: y. (x ) = Using a simple first-order difference 

i v nr i 

scheme, the equations (B.3) between points m and iih- 1 become 


.<»*» 4l . . /*IV 


<- 6 > 


and conditions (B.4) and (B.5) become 
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f t> n 1 \ 

I ) 



and 


* 


l 

I 


G 


(B.8) 


respectively. 

Our approach will be to find a method of repeatedly 
correcting a given set of trial functions until we obtain 
a set which satisfies the equation. To this end, let us 
assume a set of trial functions y^ (x) having values y^ at 
the mesh points and a trial eigenvalue X (obtained either 
from an initial guess or from the previous iteration). The 
y™ do not satisfy the equations (B.6), but we now seek an 
improved set _ 

tC- r * 

A * X 4 

Then, with the notations 
we write * 



(summation on repeated lower indices, not on repeated upper 
indices). The difference equations (B.6) are now linearized 
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by inserting the above expansions and neglecting terms of 
2 2 

order 6y i and 6X . The result, after some rearrangement, is 


(iC-ir-yiir *(«7*i=. s ‘>> s »r i 

V (B.9) 


- z^^r*'-aD-Kr v r) 

and the conditions (B.7) and (B.8) become 

s ll * (Sj-cJ-o. 


y J 


(B. 10) 


Between any two mesh points, then, we must compute the 
N x N matrices 


a ‘> ' l i ij t °>'j = [ 2 i j- ij] 


1 


and the 1 x N column vectors 


KB- 11) 


cr- kcht), er- 


With these definitions, (B.9) becomes 


+ + c t - iX + 6 t " 0. (B.12) 


Defining now the inverse a -1 . 1 ? of a ™ ; viz.. 
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(B. 13) 



and the quantities 


6 m -/ flm . /V *1 mm/m . M. ^ 

“ a ** Cj=-a ih Z ii/ E t e* , (b. i4) 


one obtains 


HT"- K + 0 * c. 


(B.15) 


We now seek a similar relation between the corrections 
at any given point and those at the first boundary point, 
thus: 


*' ~( i it <■ y'iX + £~ (B. 16 ) 

The 0. . , y™ y and e™ can be found because, clearly, 
xj x i 


4 ^ 4 

fat = 6 *i , 

v, 1 - 

c 1 

'-t , 


f. 1 


(3.17) 

and, combining (B.15) and (B.16) we 

obtain a 

set of 

recurrence 

relations : 









/m -1 





fat - 


Phi 





jrr “ 


y **”* 

•ft 

f cr 



(B.18) 

£7 - 

b7« 

M “1 

► E? 

m 




J 
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Beginning with (B.17) and applying (B.l8) repeatedly, we 
finally arrive at a relation between the corrections at 
the two ends of the interval [ 0, 1 ] : 




(B. 19) 


These constitute N linear equations relating the 2N + 1 
M+l 1 

quantities 6y_ i , dy^ Uh. In addition, we have the N + 1 
relations (B.10) and upon eliminating the 6y^ +1 from (B.19) 
with the first of (B.10), we have N + 1 equations in N + 1 
unknowns, which we then solve for the 6y^ and 6 'h. Applying 
(B.15) successively at each mesh point, we then obtain all 
the corrections 6y£ at each point. The corrections are added 
to the and A and the whole procedure is repeated. 

It may be mentioned that in difficult cases, and especially 
when the trial function is poor, it often proves expedient 
to apply not the full correction, but only a certain fraction, 
v , of it. Thus the new functions would be 


= sr * o<><i . 

After a number of iterations v may be increased, and when 
the "solution" becomes so good that the linear approximation 
(B.9) is quite good, the convergence is speeded by setting 
v = 1. In the present calculations it was found that, when 
convergence was obtained, the largest relative correction 
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to any variable at any mesh point. 



as a typical value. 

When the procedure has converged, the y^x) are periodic 

functions satisfying (B.6). If now ire repeat the process 

once more, the 6y. are just the w (x) of the Floquet theory 

i i 

(cf. Eq. 3.4), and the Floquet matrix A is evidently just 



the matrix already calculated for use in Eq. (B.19). 
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Figure 2. Time integration at R = 1000 and T = 350 showing 
the irregular modulation. 




Figure 3. Time integration at R = 10,000 and T = 3500. There 
are more oscillations in each cycle of the irregular modulation 
than in Figure 2. 
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Figure 4a. A sketch of the energy curves for B* 
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Figure 4b. A sketch of the energy curves for 



Figure 4c. A sketch of the energy curves for B* 



Figure 5 . The periods of the symmetric, asymmetric and two- 
peaked (symmetric) solutions. The period is in dynamical units 
of time. The solid curves give the periods calculated by the 
method of averaging; the dashed curves give periods at R = 100 
from the numerical relaxation solutions. The values of P ASM 
from the asymptotic theory valid as 6 -» 3/4“ are shown for 
comparison. 
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Figure 8a. Time integration of (2.8) for 6 = 0.6. The singular 
point is stable and the solution is spiralling in. 
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Figure 8b. Time integration of (2.8) for 6 = 0.625. The 
singular point is an unstable spiral and the solution enters 
a limit cycle. The limit cycle crosses OA from left to right 
leading to the sudden jump in amplitude of the periodic solution. 
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Figure 10. The quantity a - In la^ |, where is 

the Floquet multiplier of greatest modulus, for the symmetric 
asymmetric and 2-peaked (symmetric) periodic solutions at 
R = 100. Negative values of a mean that the corresponding 
periodic solution is stable; positive values, unstable. The 
dashed portion of the curve for the asymmetric solution shows 
the range in which the Floquet multipliers are complex. 





>2 


Figure 11. Several periodic solutions found at R = 100, 

6 = .748. The single-peaked solution (a) and the asymmetric 
double-peaked solution (b) are stable. The symmetric double- 
peaked solution (c) is unstable. The time axis is labeled 
in dynamical units. 




